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A formula is derived that relates the ground-state dispersion of a many-body system with the end- 
to-end distribution of paths with open boundary conditions in imaginary time. The formula does 
not involve the energy estimator. It allows direct measurement of the ground-state dispersion by 
quantum Monte Carlo methods without analytical continuation or auxiliary fitting. The formula 
is applied to the lattice polaron problem. The exact polaron spectrum and density of states are 
calculated for several models in one, two, and three dimensions. In the adiabatic regime of the 
Ifolstein model, the polaron density of states deviates spectacularly from the free-particle shape. 
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I. INTRODUCTION 

Usually, quantum Monte Carlo (QMC) methods are 
used to study ground-state, thermodynamic, or static 
properties of quantum-mechanical systems. Some im- 
portant dynamical characteristics can also be obtained 
through various forms of fluctuation-dissipation rela- 
tions. Examples of these are the superfluid fraction of 
Bose-liquids [||, the Drude weight of conductors 1^, the 
Meissner fraction of superconductors 13, and the effec- 
tive mass of defects and polarons g|| . Beyond that 
dynamical calculations are less straightforward. For in- 
stance, calculation of the excitation spectrum normally 
requires measurement of the Green's function at imagi- 
nary times and subsequent analytical continuation to real 
times. 

However, there exists one special type of excitation 
spectrum that can be measured directly by QMC. This 
is the ground-state dispersion, i.e., the total energy of the 
system Ep as a function of the total momentum P. In a 
translationally-invariant system, P is a constant of mo- 
tion, and the Hamiltonian does not mix subspaces with 
different P. Then, if a QMC is designed as to operate 
within a given P-subspace only, it may be able to ac- 
cess the ground state for the given P, thereby providing 
E-p. Not for any physical system E-p is of interest. For 
a collection of identical particles, for instance, one has 
simply E-p — P^/(2M), M being the total mass, which 
corresponds to free movement of the system as a whole. 
Positive examples include cases when the system can be 
divided into a tagged particle and an environment (usu- 
ally bosons). The best known example of this kind is 
the polaron, i.e., an electron strongly interacting with 
phonons. In this case Ep is nothing but the polaron 
spectrum. The polaron spectrum will be the main sub- 
ject of this paper. 

There exist at least two different strategies of how to 



operate within a restricted P-subspace. The first one is to 
work in momentum space and to fix the total momentum 
of the system from outset. An example of this approach 
is the diagrammatic method of Prokof 'ev and Svistunov 
[Q. In this method QMC is used to sum the entire dia- 
grammatic series for an imaginary-time Green's function 
G(P,r). Since the total momentum is an external pa- 
rameter of the series, it is possible to extract Ep from 
the T — > oo limit behavior of G, by fitting it to a single 
exponential e^^^'^ . This method is exact and universal 
but requires a separate simulation for each P-point. 

The second strategy is to work in real space but to 
use Fourier-type projection operators to project on states 
with definite P. This amounts to the free boundary con- 
ditions in imaginary-time. Usually, the projections are 
used to calculate the second derivative of the energy with 
respect to momentum (effective mass) |^,^. In Ref. |^ 
the projection was applied for the first time to the whole 
polaron spectrum. In this scheme, the ground-state dis- 
persion is measured directly, and all Ep are calculated 
simultaneously. Unfortunately, at non-zero P the weight 
of the polaron path is no longer positive definite, and one 
needs to deal with a sign-problem. It turns out, however, 
that the main idea can be reformulated in a way that 
does not require a division by the average sign, but only 
taking its logarithm. While the new formulation does not 
constitute a complete elimination of the sign-problem, it 
is more statistically stable and extends the parameter do- 
main accessible in practical simulations. Below we derive 
the new formula, discuss its properties, and apply to the 
physically interesting example of the lattice polaron. 

II. A FORMULA FOR THE GROUND-STATE 
DISPERSION 

Up to our knowledge, the projection relations required 
for our purposes were first derived by Basile p| . For 



1 



completeness a derivation is given below. Let R denote a 
many-body real-space configuration, and i? -I- r a many- 
body configuration which is a result of the parallel trans- 
port of i? by a vector r. (Note that the sum of R and 
r is only symbolic. The dimensionality of R is equal to 
the number of degrees of freedom, i.e., very large or infi- 
nite, while the dimensionality of r is the dimensionality of 
space, i.e. 1, 2, or 3.) States \R) form a complete orthog- 
onal basis, I = / dR\R){Rl and {R\R') ^ S{R - R'). 
A different basis is formed by the states |n) which are 
characterized by the definite total momentum P. One is 
interested in the projected partition function Zp which 
includes only states with the given P: 



Zp = Y,He-''"\n)Sp,p„ 

dRdR'{R'\e-^"\R) -Qp, 



(1) 



Qp = ^(i?|n)(n|i?') Sp^p„ = {R\P){P\R'), (2) 

n 

where (3 — {ksT)^^ is the inverse temperature and H is 
the full Hamiltonian. The meaning of Eq. (H) is that the 
two configurations, R and i?', have to be projected on the 
given momentum P. This is achieved as follows (below 
7i = 1 is set). Any arbitrary configuration R generates 
a set of states \Pr) = J dve-'''^''\R -f r), where V 

is the volume. Inversely, \R + r) = V'^^^ J2p e'^'"|P_R)- 
Upon projection, only P-components of both configura- 
tions survive. As a result 

gp = l I dv{R + r\P){P\R' +r) 

= T72 [drY,e'(^'-^'>{P'n\P){P\P"R') 
J p,p„ 

= 1^ j dv dr' {R + r\R' + r') e*P(>-"-') 

= iy" d(Ar)(i? + Ar|i?')e'P^'- 

= 1 y" d(Ar) e'P^"-(5((i? + Ar) - i?'), (3) 

where Ar = r — r'. Substitution in Eq. (^ and integra- 
tion over R' yields 

= y d(Ar) e*P^"-y dR{R + Ar\e-l^"\R) 

= Jd{Ar)e'^^'' JdRp{R,R + Ar;f]), (4) 

where p{R^ i?'; (3) is the full many-body density matrix. 
Next, we assume that for each P the state with the lowest 
energy Ep is non-degenerate, and in the low-temperature 



limit the projected partition function is dominated by the 
contribution from this state, Zp — > exp(— /3i?p). Now 
take the ratio of Zp and Zp=o' 



/3^oo Zp^o 



(5) 



lim 



Jd{Ai 



JPAr 



J dRp{R, R + Ar; (3) 



/d(Ar) J dRp{R, R + Ar; (3) 



where Eq is the ground-state energy. The rhs is nothing 
but the average value of cosPAr taken over the distri- 
bution p. [We have assumed that / dRp{R, R + Ar; /?) is 
an even function of Ar.] A simple formula for Ep now 
follows 



Ep - Eo = - lim iln(cosPAr) 

/3^oo (3 



(6) 



which is the main result of this section. 

Eq. (^) shows that the ground-state dispersion can 
be obtained from the end-to-end distribution of many- 
body paths. It offers a direct way of evaluating the 
ground-state dispersion by QMC methods in cases when 
p{R,R';f3) is positive-definite. However, the QMC pro- 
cess must be organized in a special way, as apparent from 
Eq. (|^). It must generate only such paths whose end con- 
figurations at imaginary time t = (3 are exact images of 
the end configurations at r = except for a parallel trans- 
port by an arbitrary vector Ar. It is allowed to change 
Ar, make simultaneous changes of both end configura- 
tions, and to make arbitrary changes of paths at internal 
times < r < /3, but the end configurations must always 
be kept identical up to a shift. It is important that this 
restriction affects neither the ergodicity nor the applica- 
bility of the Metropolis algorithm. 

Formula (^ involves only one measures quantity, 
(cosPAr), instead of two in the previous formulation. 
Moreover, it does not require division by the measured 
quantity, but only taking its logarithm. Additionally, 
Eq. (H) provides the difference between two large num- 
bers, Ep and i^o, and a large cancellation of errors may 
occur. This makes Eq. (^) much more stable statistically 
than explicit energy estimators. 

On the other hand, at small temperatures, the aver- 
age cosine becomes exponentially small and it cannot be 
measured reliably. This reflects the fact that configura- 
tions with Ep — Eq ^ ksT are very rare because of the 
Boltzmann factor. Thus, the present method is limited 
to excitation energies of the order of several fc^T. 



III. APPLICATION TO THE LATTICE POLARON 

We now demonstrate the practical importance of 
Eq. (^ on the model problem of lattice polaron, which 
is often considered as a paradigmatic example of a parti- 
cle strongly interacting with a boson field. We consider 
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a hypercubic lattice with the nearest-neighbor hopping, 
dispersionless phonons, and the "density-displacement" 
electron-phonon interaction. The model Hamiltonian 
reads 



H = - 



(„„') 



E 



huj ^ bl^bm. 



(7) 



Here t is the hopping amplitude (it will be used as the 
energy unit) , uj is the phonon (oscillator) frequency, is 
the internal coordinate of the mth oscillator, and /m(n) 
is the force between mth oscillator and the particle at 
site n (/ is a function of distance |m — n| only). The 
model is parametrized by the dimensionless frequency 
uj — huj/t and by the dimensionless coupling constant 
A = Em /m(0)]/(2Mw^ D), where M is the mass of the 
oscillator and D is the half-bandwidth of the bare band. 
(For an isotropic band with nearest-neighbor hopping, 
D = zt, z being the number of neighbors.) 

For the polaron problem, a many-body configuration 
R is specified by the position of the electron r and os- 
cillator displacements ^m. Making use of the Feynman's 
idea of analytic integration over jsj the problem is 
reduced to a single-particle system with retarded self- 
interaction. The latter can be simulated exactly, using 
the continuous-time representation of polaron paths 
The resulting algorithm Q is very efficient and allow 
accurate determination of the ground-state energy and 
effective mass of the polaron for a wide class of mod- 
els. It this paper, it will be shown that the method also 
produces accurate polaron spectra, when combined with 
Eq. (|^). There are other reasons why the polaron is an 
ideal system to try formula (^. First, due to a constant 
phonon frequency, excited states are, at any P, separated 
from the restricted ground state by a finite energy gap 
huj. Therefore, instead of performing numerically the 
limit procedure to /3 = cx3, one can study the system at 
finite f3, provided exp(/3?iw) 3> 1 and the contribution 
from excited states is negligible. Second, by increasing 
the coupling constant A one can always decrease Ep — Eq, 
i.e., substantially increase (cosPAr), and stabilise the 
sumulations. Third, the polaron momentum P is not 
a parameter of simulations. This implies that statistics 
can be collected for all momenta simultaneously. In other 
words, the whole polaron spectrum is measured in a sin- 
gle QMC run. This will enable us to calculate for the 
first time exact polaron densities of states. 

We begin with the simplest Holstcin model which has 
local electron-phonon interaction, /m(n) — nSmii- In one 
dimension, the polaron spectrum has been extensively 
studied by exact diagonalisation (nOUO], strong-coupling 
perturbation jl^ , and variational [|13[ techniques. Our 
QMC data for the one-dimensional Holstein model are 
shown in Fig. ^. The most interesting feature of the 
spectrum is its non-cosine shape in the adiabatic regime 
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FIG. 1. Polaron spectrum in the one-dimensional Hol- 
stein model, normalized to the bandwidth W = Et^ — Eq. 
Triangles: Cj = 1.0, A — 2.0 [for these parameters the polaron 
ground-state energy Eo = —4.38(1) (in units of t), bandwidth 
W = 0.1243(2) (in units oit), effective mass m* = 10.0(1) (in 
units of mo — fi^ /{2ta'^)).] Diamonds: lj — 1.0, A — 2.5 
[Eo = -5.26(1), W = 0.0437(3), m* = 34.5(3)] . Cir- 
cles: Lj = 10.0, A = 10.0 [Eo = -20.35(1), W = 0.543(2), 
m* = 6.06(2)] . Squares: lo = 10.0, A = 20.0 [Eo = -40.08(1), 
W = 0.0739(2), m* = 47.6(1)] . 



LU < 1.0 (triangles and diamonds). At large momenta, the 
spectrum is more flat than at small momenta. The nature 
of this flattening was understood a long time ago . In 
the weak-coupling limit, the free-particle state hybridizes 
with the single-phonon state and creates a mixed ground 
state which is free-particle-like at small P and phonon- 
like at large P, hence the weak dispersion. With increas- 
ing A, the free-particle state is replaced with a polaron 
state with an effective mass m*, which now hybridizes 
with the single-phonon state, still leading to a more flat 
dispersion at large momenta. The flattening effect weak- 
ens with growing A and (D because both processes increase 
the energy separation of the two hybridizing states. In 
recent years the flattening of the polaron spectrum was 
observed in numerical studies 10 1^. 

Our Quantum Monte Carlo data fully confirm the pre- 
vious analytical and numerical results, see Fig. 0. We 
found that the spectrum shape is more sensitive to the 
phonon frequency than to the coupling constant. At a 
small frequency <D — 1.0, the increase of the coupling 
constant from A = 2.0 to A = 2.5 results in a 3.5-times 
increase of the effective mass, and in a 2.8-times drop 
of the bandwidth, yet the spectrum shape changes only 
slightly, cf. triangles and diamonds in Fig. |l[ At the 
same time, a simultaneous 10-times increase of A and uj 
results in a similar, 4.8 times, increase of effective mass 
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but brings the spectrum shape very close to cosine, cf. 
triangles and squares in Fig. |^. Again, a doubling of A 
strongly affects Eq and W but not the spectrum shape, 
cf. circles and squares in Fig. |^. 
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FIG. 2. Spectrum of the two-dimensional Holstein model 
in the adiabatic regime, uj = 1.0, A = 1.4. [Eq = —6.12(3), 
m* = 8.7(1).] 

It is instructive to compare the exact QMC results with 
the Lang-Firsov (LF) approximation which is be- 
lieved to be the correct description of the polaron in the 
antiadiabatic regime (high phonon frequency). The LF 
formula for the spectrum reads 



Ep- Eo = 2ie-^^/^(l - cosP) 



(8) 



which also implies the relation between the renormalized 
mass and bandwidth 



mo 



Wo 
W 



(9) 
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where — 2zt is the bare bandwidth and 
Ti^/{2ta?) is the bare mass, a being the lattice con- 
stant. For uo = 10.0 and A = 10.0 QMC results are 
W = 0.543(2) t and to* = 6.06(2) too while LF yields 
Wlf = 0.541 1 and m^^ = 7.39 too. For u = 10.0 and 
A = 20.0 one has W = 0.0739(2) t and to* = 47.6(1) too 
from QMC and Wlf = 0.0733 1 and mlp = 54.6 too from 
LF. One can see that LF predicts very accurate values of 
the polaron bandwidth. This fact was established in the 
previous studies of the Holstein model ||l^,|l^. On the 
other hand, LF slightly overestimates the polaron mass. 
This is due to small deviations from the cosine shape, 
still present in the true spectrum at these model param- 
eters. Still, the LF masses are reasonably close to the 
exact ones, and the agreement improves with the further 
increase of u and A. 

Consider now the two-dimensional Holstein model. 
The only exact polaron spectra published so far were 



calculated by Wellein, Fehske, and Loos with the exact 
diagonalization method . These authors found a flat- 
tening of the spectrum in the outer part of the Brillouin 
zone, even stronger than in the one-dimensional case. We 
checked that for the model parameters used in [|ll|, for- 
mula (^ yields precisely the same values of E^ — Eq as 
the exact diagonalization method. A definite advantage 
of the present method is that it allows simultaneous cal- 
culations at any desired number of P-points, while exact 
diagonalization studies are limited to a small number of 
P-points due to the finite size of the clusters. On the 
other hand, the QMC method is limited to the condition 
W <^ huj, which prevents us from studying the weak- 
coupling regime and such an interesting phenomenon as 
the limit point of the polaron spectrum ||l4[] (The latter 
is possible with the diagrammatic QMC0.) Figure ^ 
shows our QMC data for a new set of parameters in the 
adiabatic regime, w = 1.0, A = 1.4, where 30 P-points 
have been used to represent the spectrum. One can see 
that the dispersion is indeed weak for |P| > 7r/2, i.e. 
in the larger part of the Brillouin zone. Note, that at 
these parameters the polaron bandwidth is reduced by 
8t/0.12i = 67 times, and the mass enhancement is 8.7, 
so we are already in the small polaron regime. Yet, the 
spectrum shape is profoundly non-cosine. With increas- 
ing A, it will be approaching the cosine shape, but this 
is expected to happen only at such large A where the 
polaron is very heavy and is easily localized. 
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FIG. 3. Density of states of the two-dimensional Hol- 
stein model in the adiabatic regime, to = 1.0, A = 1.4. 
[£o = -6.12(3), m* = 8.7(1).] 

In two dimensions, the volume of the outer part of 
the Brillouin zone is larger than that of the inner one. 
Therefore the linear representation of the spectrum, like 
in Fig. 1^, does not fully convey the changes in the band 
structure, caused by the flattening of the spectrum. The 
proper physical quantity which takes into account all the 
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states of the Brillouin zone is the density of states (DOS). 
The QMC method, coupled with formula provides 
the unique opportunity to calculate polaron DOS exactly, 
since it allows a simultaneous measurement of the whole 
spectrum. In this work, the two-dimensional Brillouin 
zone was divided in 200^ points at which the spectrum 
was measured. In the end, the total of 40 000 polaron en- 
ergies were distributed over 50 energy intervals between 
and W. The resulting DOS for Q = 1.0, A = 1.4 (the 
same parameters as in Fig. ^ is shown in Fig. ^. One 
can see that the effect of the spectrum flattening is indeed 
quite dramatic. The upper half-band is jammed into a 
narrow, 0.015 f-width, energy interval, thereby increasing 
DOS at the top of the band to « 50 times the DOS at the 
bottom of the band. The van Hove singularity is shifted 
from the middle to the top of the band. The lower half of 
the band contains only 13% of all states. Overall, DOS 
looks qualitatively different from the free-particle one. 
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FIG. 4. Spectrum of the two-dimensional Holstein 
model in the anti-adiabatic regime. Oj = 8.0, A = 8.0. 
[Eo = -32.16(1), m* = 38.4(1).] 

As in the one-dimensional case, the band structure ap- 
proaches the free-particle-like as the phonon frequency 
increases. As an example, we considered a frequency 
equal to the bare bandwidth, uj = 8.0, and A = 8.0. 
(This value of the coupling constant was chosen to have 
the polaron bandwidth, VF, close to the previously con- 
sidered adiabatic case.) The polaron spectrum and DOS 
are shown in Figs. ^ and ^, respectively. Deviations from 
the free-particle behavior are still visible, but they are 
small and not qualitative. DOS at the top of the band 
is just w 1.5 times larger than at the bottom, and the 
singularity appears close to the middle of the band. The 
"normalization" of the spectrum at large Co is quite un- 
derstandable. With increasing phonon frequency, the re- 
tardation effects become less important, the lattice defor- 
mation more readily follows the particle movement, and 
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FIG. 5. Density of states of the two-dimensional Hol- 
stein model in the anti-adiabatic regime, u) = 8.0, A = 8.0. 
[Eq = -32.16(1), m* = 38.4(1).] SmaU fluctuations of DOS 
are artifacts of a finite-mesh integration over the Brillouin 
zone and not statistical errors. 

the whole complex behaves more like a free particle with a 
renormalized hopping integral. The Lang-Firsov formula 
(|) predicts Wlf = 0.1465 i and m\p = 54.6 mo which 
is to be compared with the QMC results W = 0.1510(3) t 
and m* = 38.4(1) toq. Again, the LF approximation 
yields the correct bandwidth but overestimates the ef- 
fective mass by some 40%. 

The most spectacular transformation of DOS occurs 
in the three-dimensional Holstein model. In three di- 
mensions, the volume of the outer part of the Brillouin 
zone is much larger than that of the inner part, and it 
should completely dominate the total DOS. We do not 
show the polaron spectrum which is not very informa- 
tive. Density of states was calculated by measuring the 
polaron spectrum at 60^ points of the full Brillouin zone, 
and then distributing them among 50 energy intervals 
between and W. DOS in the adiabatic regime, cD = 1.0 
and A = 1.2, is shown in Fig. ||. The states of the flat part 
of the spectrum form a massive peak at the top of the 
band. The width of the peak is about 10% of the total 
bandwidth. DOS at the bottom of the band is negligible, 
the capacity of the lower half-band is less than 1% of the 
total number of states. The two van Hove singularities 
are not visible, at least on the chosen level of energy res- 
olution, they are absorbed into the peak. Overall, DOS 
is completely different from the free-particle one. Should 
the three-dimensional Holstein model with such parame- 
ters exist in nature, an extreme care would be necessary 
in interpreting experimental data. In any real material, 
the lowest states would likely be localized, and any re- 
sponse to an external perturbation would be dominated 
by the peak. Then, for instance, fltting to a free-particle- 
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FIG. 6. Density of states of the three-dimensional Hoi- 
stein model in the adiabatic regime, u) — 1.0, A = 1.2. 
[Eo = -7.75(4), m* = 6.2(2).] 

like form of DOS would lead to wrong estimates of the 
coupling constant and other errors. 

As before, the band structure returns to the free- 
particle shape in the anti-adiabatic regime. We consid- 
ered the case of the phonon frequency being equal to 
the bare bandwidth, uj = 12.0, and coupling constant 
A — 10.0, when the polaron bandwidth is close to the 
just considered adiabatic case, see Fig. |^. Although still 
distorted, the DOS shape is close to the free-particle one, 
with the square-root behavior at the top and the bottom 
of the band, and with two van Hove singularities fully 
developed at the "right" places. The polaron bandwidth 
is W = 0.0827(2) t, which is in good agreement with the 
Lang-Firsov value Wlf = 0.0809 1, while the polaron 
mass m* — 112(1) toq is 24% lighter than the LF mass 
m}^p = 148 toq. 

The polaron QMC algorithm of Ref. |^ is not limited to 
the Holstein model. In fact, it allows studies of arbitrary 
forms of the electron-phonon interactions (of the density- 
displacement type), and arbitrary forms of the particle 
kinetic energy. Combined with Eq. (^, it provides an 
efficient and exact way of calculating the band structure 
of a whole class of polaron models. As possibilities are 
numerous, we have chosen to illustrate the point on two 
particular examples. 

The first example is the anisotropic two-dimensional 
Holstein model with (D = 1.0, A — 1.4 (these param- 
eters are the same as in Figs. ^ and ||), and the bare 
anisotropy ratio ty/t^ = 0.2. For a free particle with such 
an anisotropy, the saddle points at (±7r,0) and (0, ±7r) 
have different energies, which results in two singularities 
in DOS, positioned symmetrically with respect to the 
center and edges of the band. Polaron DOS, calculated 
by QMC, is shown in Fig. ||. The flattening effect cre- 
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FIG. 7. Density of states of the three-dimensional Holstein 
model in the anti-adiabatic regime, lj = 12.0, A = 10.0. 
[Eo = -60.12(2), m* = 112(1).] Small fluctuations of DOS 
are artifacts of a finite-mesh integration over the Brillouin 
zone and not statistical errors. 

ates a strong peak at the top of the band which absorbs 
the higher-energy singularity [at (±7r, 0)]. At the same 
time, the second singularity is still clearly visible. Now 
it appears in the vicinity of the middle of the band. 

The second example is the two-dimensional polaron 
model with long-range electron-phonon interaction [com- 
bine with the Hamiltonian (^]: 

7\ I 1\3/2 ' (l*^) 

(|m — n|"^ -I- ly/-' 

where the distance |m — n| is measured in lat- 
tice constants. For this form of the force, A — 
1.7A2k'^ /{2Mlj'^D). The model describes a two- 
dimensional particle interacting with a parallel plane of 
ions vibrating perpendicular to the plane. It was pro- 
posed in Ref. [^ , where it was used to model the interac- 
tion of holes doped in copper-oxygen planes, with apical 
oxygens in the layered cuprates. It was found that this 
long-range (Frohlich) polaron is much lighter than the 
short-range Holstein polaron. Here we present the den- 
sity of states for uj — 1.0 and A = 2.75, see Fig. |9[ DOS 
shape is close to the free-particle one, with a single, well- 
developed singularity in the middle of the band. Note, 
that we are in the adiabatic regime, at the same frequency 
w = 1.0 where the two-dimensional Holstein polaron has 
a very distorted DOS, cf. Fig. ^. The comparison of 
Figs. 1^ and ^ shows that a long-range electron-phonon 
interaction plays the same role as increasing phonon fre- 
quency, as far as the flattening effect is concerned. 
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FIG. 8. Density of states of the two-dimensional 
anisotropic Holstein model. iJo — 1.0, A = 1.4, ty/tx = 0.2. 
[Eo = -3.987(3), ml = 3.44(1) mo^, m% = 17.54(3) mo^.] 

IV. DISCUSSION AND CONCLUSIONS 

The main message of this paper is that a path-integral 
imaginary-time Quantum Monte Carlo is quite capable 
of direct measuring of real-time spectra. By relaxing the 
boundary conditions in imaginary time, one allows many- 
body paths to have arbitrary real-space shift Ar. Then 
the Fourier transform projects out configurations with a 
certain total momentum P. The result, expressed by the 
formula (^, is the ground-state energy as a function of 
the total momentum. The physical system, where such a 
ground-state dispersion is of interest, should be carefully 
chosen. 

The role of temperature in this process is two-fold. On 
one hand, temperature should be made as low as possi- 
ble to exclude the contribution from the excited states 
with the same P. On the other hand, only states with 
E-p — Eq ~ ksT are excited in the system. This means 
that the corresponding configurations will be generated 
by the QMC process in amounts sufficient for a good 
statistics. For higher-energy states, configurations will 
be exponentially rare. This is the reason for the average 
cosine in Eq. (^) to become exponentially small in the 
low-temperature limit. In this case, the measurement 
process will be statistically unstable. Thus, the temper- 
ature should be of the order of the energy interval one 
is interested in. The two conditions on the temperature 
can be reconciled if the energy scale of the ground-state 
dispersion is much smaller than the energy gap between 
the ground state and the excited states within the same 
P-sector. This is realized in the polaron system where 
one has W -^Tiu! for a wide range of parameters. If this 
condition is not satisfied, Eq. (o) will be measuring the 



FIG. 9. Density of states of the two-dimensional model 
with long-range electron-phonon interaction. lj = 1.0, 
A = 2.75. [£^0 = -11.83(1), m* = 20.1(1).] Small fluctu- 
ations of DOS are artifacts of a finite-mesh integration over 
the Brillouin zone and not statistical errors. 

difference of projected free energies rather than that of 
ground-state energies. 

Eq. ^ shows that in a many-body system there exists 
a general and simple relation between the ground-state 
dispersion and the end-to-end distribution of imaginary- 
time paths. It does not involve the energy estimator, 
although the latter is required for the separate calcula- 
tion of Eq and Ep. This might be useful in cases when 
the evaluation of the energy estimator is computationally 
costly. There are other computational advantages. First, 
Eq. (^) involves only the logarithm of one measured quan- 
tity, the average cosine. Second, Eq. (^) calculates the 
difference of two energies both of which may be large. In 
the polaron problem, typical energies are of the order of 
a few t but the bandwidth isW^ 0.1 1. Both Eq and Ep 
can be calculated with typical accuracy 0.3 — 0.5 %. This 
may result in a sizable error in their difference if the two 
energies are calculated separately and then subtracted. 
Eq. (|^) produces much more stable energy differences be- 
cause of the large cancellation of errors between Ep and 
Eq. In all the spectra presented in this paper, Figs. |^, 
^, |[ the statistical errors are smaller that symbols rep- 
resenting the data. Finally, the whole dispersion (as well 
as Eq, and derivatives at any P-point) can be calculated 
during a single QMC run. This property of Eq. (^) also 
allows fast computation of the density of states. 

To demonstrate the practical usefulness of Eq. (|^), 
we have combined it with an exact continuous-time al- 
gorithm for the lattice polaron and calculated first 
detailed polaron spectra and densities of states in two 
and three dimensions. Although our method of calculat- 
ing the spectrum is limited to the condition W <^ fiuj, 
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i.e., to the intermediate and strong coupling, this is 
the most physically interesting regime. Together with 
the weak- and strong-coupling perturbation expansions, 
the exact diagonalization, density-matrix renormaliza- 
tion group variational, and diagrammatic QMC 

techniques, the present method covers the whole parame- 
ter range of the polaron problems. One can say, that the 
problem of calculating the exact polaron spectrum has 
found its solution. 

Apart from being a test for Eq. (^) , the lattice polaron 
still has a considerable interest on its own. We have seen 
in this paper that the polaron spectrum in the adiabatic 
limit of the Holstein model is generically non-cosine, as 
was predicted theoretically and recently observed numer- 
ically. The flattening has been found to continue well into 
the strong-coupling regime, where the polaron mass is a 
few dozens and approaches a hundred. The same conclu- 
sion was reached previously in [pT| . For physical applica- 
tions the most interesting regime is when polarons arc not 
very heavy and can be mobile. We conclude that the non- 
cosine spectrum is typical for the Holstein model in the 
physically relevant parameter region, i.e., small phonon 
frequencies and intermediate couplings. 

A surprise finding has been the great extent at which 
the flattening changes the band structure in high dimen- 
sions. Density of states is completely changed, van Hove 
singularities are shifted, low-lying states are almost ir- 
relevant, relations between the effective mass and the 
bandwidth is broken. We have seen that the combination 
of dimensionality, phonon frequency, coupling strength, 
and anisotropy may produce densities of states of various 
shapes. In such a situation, one should be careful about 
the interpretation of any experimental data, like the esti- 
mation of A or polaron mass from the bandwidth or from 
the location of singularities. The use of a "naive" model 
band structure would lead to wrong conclusions. 

We have also found that in all dimensions the po- 
laron band structure becomes free-particlc-like with in- 
creasing phonon frequency In particular, the spectrum 
approaches the cosine shape and the effective mass ap- 
proaches the inverse of the bandwidth. Moreover, numer- 
ical values of W are well described by the Lang-Firsov 
approximation. Thus, our QMC data support the LF 
approximation as the right description of the polaron in 
the antiadiabatic regime. At the same time, we have 
found that LF could still overestimate the polaron mass 
by a few dozen percent in cases where the bandwidth is 
predicted correctly. 

Finally, we considered a long-range electron-phonon in- 
teraction and found a free-particle-like band structure 
even in the adiabatic regime. In Rcf. it was found 
that in the adiabatic regime polaron masses arc well de- 
scribed by the LF approximation. Two conclusions fol- 
low from these facts. First, all the unusual properties 
of the Holstein model caused by the flattening, may be 
specific to the local electron-phonon interaction and may 



not be generic polaron properties. Second, the long-range 
electron-phonon interaction on a lattice seems to have 
the same effect on the band structure as the increasing 
phonon frequency. Although it is clear intuitively that 
a long-range interaction leads to higher mobility of the 
lattice deformation, details of this mechanism are yet to 
be fully understood. 
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